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We discuss dynamics of approximate adiabatic invariants in several nonlinear models being related 
to physics of Bose-Einstein condensates (BEC). We show that nonadiabatic dynamics in Feshbach 
resonance passage, nonlinear Landau-Zener (NLZ) tunnelling, and BEC tunnelling oscillations in a 
double- well can be considered within a unifying approach based on the theory of separatrix crossings. 
The separatrix crossing theory was applied previously to some problems of classical mechanics, 
plasma physics and hydrodynamics, but has not been used in the rapidly growing BEC-related field 
yet. We derive explicit formulas for the change in the action in several models. Extensive numerical 
calculations support the theory and demonstrate its universal character. We also discovered a 
qualitatively new nonlinear phenomenon in a NLZ model which we propose to call separated adiabatic 
tunnelling (AT). 



I. INTRODUCTION 

Adiabatic invariance [l[ is very important in many 
fields of physics. In the last decade, there has been a 
great deal of interest in physics of Bosc-Einstcin con- 
densates @, H, H], H, H, 0] among scientists from several 
scientific fields. Presently BEC research is at the cross- 
ing point of AMO science, statistical mechanics and con- 
densed matter physics, nonlinear dynamics and chaos. 
The discussion we present here is related to interplay 
between nonlinearity and nonadiabaticity in BEC sys- 
tems. The relation between quantum transitions and 
change in the classical action of a harmonic oscillator 
has long been known [H, [j| . BEC bring nonlinearity into 
a quantum world. BEC dynamics can often be described 
within the mean-field approximation; finite-mode expan- 
sions produce nonlinear models where a variety of phe- 
nomena common to classical nonlinear systems happen. 
We consider two kinds of nonlinear phenomena here: de- 
struction of adiabatic invariance at separatrix crossings 
and probabilistic captures in different domains of phase 
space. 

A conceptual phenomenon of classical adiabatic the- 
ory is destruction of adiabatic invariance at separatrix 
crossings which is encountered, in particular, in plasma 
physics and hydrodynamics, classical and celestial me- 
chanics [IS m, E3, M, M M M Ei M M MMM- 

The phenomenon is very important for BEC physics: we 
consider here nonlinear two-mode models related to tun- 
nelling between coupled BEC in a double well [23| , non- 
linear Landau-Zener tunnelling 0, |25| . Feshbach res- 
onance passage in atom- molecule systems [III H3, [Hj]. 
Nonlinear two-mode models were extensively studied pre- 
viously (sometimes beyond the mean-field approxima- 
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^nj23|j24U25 j29U3fj, |31|, M \M M \M M \M M 
HI SO, EH El El El), and destruction of adiabatic- 
ity was discussed already in [2I [H |2(| , still there are 
regimes of motion that were not analyzed in these pa- 
pers from the point of view of nonadiabatic behaviour, 
that is, when initial populations of both modes are not 
zero (or very small), but finite. We presented some of 
our results on that theme in [H [Hj]; nevertheless de- 
struction of adiabatic invariance has not been studied 
systematically in BEC-related models yet. Action is an 
approximate adiabatic invariant in a classical Hamilto- 
nian system that depends on a slowly varying parameter 
provided a phase trajectory stays away from separatrices 
of the unperturbed (frozen at a certain parameter value) 
system. If this condition is not met, adiabaticity may 
be destroyed. As the parameter varies, the separatrices 
slowly evolve on the phase portrait. A phase trajectory 
of the exact system may come close to the separatrix 
and cross it. Thegeneral theory of the adiabatic sep- 
aratrix crossings [lfj predicts universal behavior of the 
classical action (described in a greater detail in the main 
text). In particular, at the separatrix crossing the ac- 
tion undergoes a quasi-random dynamical jump, which 
is very sensitive to initial conditions and depends on the 
rate of change of the parameter. The asymptotic for- 
mula for this jump was obtained in [ill . |12| . 1 1 31 ] - Later, 
the general theory of adiabatic separatrix crossings was 
also developed for slow- fast Hamiltonian systems [HEt]], 
and was applied to certain physical problems (see, for 
example, [la. [l8l [l^. I20I]). It was also noticed that non- 
linear Landau-Zener (NLZ) tunnelling models constitute 
a particular case for which the general theory can be 
applied [H[28[. Beside the quasi-random jumps of adia- 
batic invariants, there is another important mechanism of 
stochastization in BEC-related models: scattering on an 
unstable fixed point with a capture into different regions 
of phase space after a separatrix crossing [H El, E3- 
Here stochastization happens due to quasi-random split- 
ting of phase flow in different regions of phase space at 



2 



the crossing. Rigorous definition of such probabilistic 
phenomena in dynamical systems was given in [48j . The 
probabilistic capture is important in problems of celes- 
tial mechanics [lOj, but it was also investigated in some 
problems of plasma physics and hydrodynamics [l6| , op- 
tics [Hj], classical billiards with slowly changing param- 
eters and other classical models As shown in (47| . 
the combination of the two phenomena leads to dephas- 
ing in dynamics of globally coupled oscillators modelling 
coupled Josephson junctions. However, it seems that the 
probabilistic capture mechanism was not discussed at all 
in relation to BEC models yet. We discovered that in a 
nonlinear Landau-Zener model such mechanism may take 
place, and it leads to a new phenomenon (in the context 
of the model) that we propose to call separated adiabatic 
tunnelling. 

Let us review the models being considered in the 
present paper in more concrete terms. The nonlinear 
two-mode model introduced in [23[ describes BEC tun- 
nelling oscillations in a double- well as that of a non-rigid 
pendulum. In the case of an asymmetric double-well, the 
effective classical Hamiltonian is: 

H = -Sw+^ \/l - w 2 cos 6, (1) 

where w, 9 are the population imbalance and phase differ- 
ence between the two modes, parameters S and A repre- 
sent potential difference between the wells and nonlinear- 
ity, correspondingly. The same Hamiltonian appears in 
a nonlinear Landau-Zener model [U . As one slowly 
sweeps 6, say, from a large positive to a large negative 
value, a change in the mode populations is determined by 
the change in the classical action (since at large |<5| classi- 
cal action depends linearly on w) . This provides interest- 
ing link between fundamental issue of classical mechanics, 
dynamics of approximate adiabatic invariants (classical 
actions) , and nonadiabatic transitions in quantum many- 
body systems. The dynamics of classical actions in non- 
linear systems is, however, a very complicated issue [Icj . 
Some analysis of the NLZ model was done in [3, [H[ . In 
[25l | so-called subcritical (A < 1), critical (A = 1), and 
supercritical (A > 1) cases were defined. However, only 
the case of zero initial action was considered, that is a 
vanishingly small initial population in one of the states. 
We concentrate on the case of finite initial action, and su- 
percritical case. In the supercritical case, the most strik- 
ing phenomenon from the point of view of physics is the 
so-called nonzero adiabatic tunnelling (nonzero AT). In 
terms of the theory of separatrix crossings, it is caused by 
the geometric jump in the action at the separatrix cross- 
ing. Mathematically, it is a very simple issue: as a phase 
point leaves a domain bounded by a separatrix of the un- 
perturbed system and enters another domain, its action 
undergoes a " geometric" change proportional to the dif- 
ference in areas of the two domains [491 ]. The gist of sep- 
arated adiabatic tunnelling (seoarated AT) that we found 
is as follows. The separatrix divides the phase portrait 
on 3 domains Cn. 2,3 (Fig. (HI). In case at the moment 
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leaving the third domain G3 (with decreasing area) can 

be captured in either of the two growing domains. Bunch 

of trajectories with close initial actions h will be "split- 
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ted" on two bunches with two different final actions I f . 
It is possible to calculate probability for a phase point 
to come to either of the two bunches (we calculated it 
in Appendix B and compared analytical prediction with 
numerical result in Fig. 11). Possible physical applica- 
tions of nonzero AT phenomenon has been extensively 
discussed (for example, [24], [H| : wavepacket in an accel- 
erated optical lattice should undergo nonzero tunnelling 
in adiabatic limit when nonlincarity is large enough, al- 
though no experimental evidence is available yet). In 
relation to BEC oscillations in asymmetric double-well, 
corresponding physical effect is (obvious) drastic change 
in the amplitude of oscillations when regime of motion is 
changed from self-trapped to complete oscillations due to 
slow change of parameters. In the case of separated AT, 
the effect for the asymmetric well may look like this: the 
asymmetry between the wells is slowly changed; regime 
of motion is changed from self-trapped to complete os- 
cillations and then back to self-trapped. But final state 
is " splitted" : a system has " choice" of two different fi- 
nal states. For a set of experimental realizations with 
close initial conditions, one can define a "probability" 
for a system to come to either of the two states. While 
such experimental realization seems to be even less real- 
istic than nonzero AT, conceptually it is a very interest- 
ing phenomenon worth discussing: the "probability" is 
of purely classical origin. Analogous interpretation can 
be done for BEC experiencing NLZ tunnelling in optical 
lattice. Although the phenomenon looks similar to the 
nonzero AT described in [24], [25| , its mathematical back- 
ground is very much different and not so straightforward; 
it is a particular case of probabilistic phenomena in dy- 
namical systems defined in [48| . 

We also derive a formula for the jump of the adia- 
batic invariant (Eq. (|2"3"|) ) in the symmetric well case 
(5 = 0) and check it numerically (for the asymmetric 
case, the corresponding formula has both terms of order 
e and elne). As physical application of this jump, one 
can imagine an experiment with BEC oscillations in a 
double-well, with the potential barrier between the wells 
being slowly raised and then slowly decreased back to 
its initial position. The system will not return back to 
its initial state. Within the mean-field two-state model, 
the difference between the initial and final oscillations 
is caused by the change in the adiabatic invariant (of 
course, in a real system many other complications arise). 
Such kind of experiments are feasible [50, [iLfl] ■ 

Similar nonadiabatic phenomena arise in coupled 
atom-molecular systems. Here, in the mean-field limit 
it is possible to construct two-mode models based on the 
all-atom and all-molecule modes, and their coherent su- 
perpositions. The two-mode model describing a degener- 
ate gas of fermionic atoms coupled to bosonic molecules 
was considered in [2(| H3, [28| (the same model enables 



3 



(a) 0>) (c) (&) 




t 2k jc 2* n 2x it 2* 



FIG. 1: Phase portraits of the Hamiltonian (gj with A = 0. From left to right: S = 10, y/2, 1, 0, -1, — \/2, -5, -50. Stars (bold 
dots): unstable (stable) fixed points. See detailed discussion in [27l. |28| 



to describe coupled atomic and molecular BECs, so we 
call it 2-mode AMBEC model). The system is reduced 
to the classical Hamiltonian 

H = -8(t)w + (1 - w) VTT^cos 6, (2) 

where w denote population imbalance between atomic 
and molecular modes, and S is the (slowly changing) de- 
tuning from the Feshbach resonance. As 5 sweeps from 
large positive to negative values, the system is transferred 
from the all-atom w = 1 mode to the all-molecule w = — 1 
mode. The final state of the system contains the non- 
zero remnant fraction, which can be calculated as change 
in the classical action in the model (PJ), and scales as a 
power-law of the sweeping rate. The model was intro- 
duced in [2(| in the attempt to describe recent exper- 
iments on Feshbach resonance passage (EH, [H, [H, \5l\ , 
and some power laws were calculated there and compared 
with experimental data. For the case of nonzero initial 
molecular fraction, the power-law was also calculated in 
[27l |28| according to the general theory. We carefully 
check numerically this (linear) power law in Section IIB. 
It is important to note that the model give 100% con- 
version efficiency in the adiabatic limit, while in the ex- 
periments finite conversion efficiency has been seen. In 
Section IIC we present brief analysis of a more general 
model, which have an analog of nonzero AT (leading to 
finite conversion efficiency in the adiabatic limit). In the 
more general version, s-wave interactions were taken into 
account, so the Hamiltonian looks like 

H = -Sw + Xw 2 + (1 - w)y/l + w cos 6, (3) 

Here, the phase portraits can have more complicated 
structure, and the passage through the separatrix can be 



accompanied by the geometric jump in the action, lead- 
ing to a non-zero remnant fraction even in the adiabatic 
limit. 

In Section III, the nonlinear two-mode model JT]) for 
two coupled BECs is considered. For brevity, we call 
this model 2-mode atomic BEC (ABEC) model. The 
separated AT is demonstrated in the end of the Section. 

Section IV contains concluding remarks. In the Ap- 
pendix A we described the adiabatic and improved adia- 
batic approximations. In the Appendix B, derivation of 
formula for probabilities in separated AT is presented. 

The most interesting new results of the paper are: 

1) Extensive numerical tests of the formula (fT2"|) for the 
dynamical jump in the adiabatic invariant of the atom- 
molecule system (Section lib). The formula is based on 
Eq.tflO]) being obtained elsewhere [27|,[28j]. 

2) Suggested mechanism (analog of nonzero AT, Sec- 
tion lie, Fig. [J) leading to finite conversion efficiency in 
atom-molecule systems due to a geometric jump in the 
action. 

3) Analytical derivation of the explicit expression (|2"3"|) 
for the jump of adiabatic invariant of the symmet- 
ric ABEC model and its numerical test (Section Illb, 
FiglinD. 

4) Discovery of new phenomenon in nonlinear Landau- 
Zener model: separated AT. Analytical calculation of 
probabilities related to this tunnelling (B32-34) and its 
numerical test (Section IIIc, Fig. ITTj) . 

The main result is demonstration of usefulness of sepa- 
ratrix crossing theory in a variety of BEC-related models. 

In order to keep the paper compact, we do not present 
here comparison with quantum many-body calculations, 
but consider only mean-field models. 
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FIG. 2: Phase portraits of the Hamiltonian (j4| with A < (A = —0.5). From upper left to bottom right (a-1): 
5 = 5.0, 1.0, 0.53, 0.5, 0.45, 0.44, 0.4, 0, —0.5, —2.2, —5, —50. In (c)-(f), the separatrix divide phase portraits on three domains 
Gi,2,3 (G2 is adjacent to the segment w = 1, Gi is adjacent to w = —1, G3 is the loop in between). Starting with small initial 
action at w ~ 1, a phase point undergoes a geometric jump in the action in addition to a dynamical jump. This leads to analog 
of nonzero AT and finite conversion efficiency in the adiabatic limit. 
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FIG. 3: Graphical solution of the equation ©. The line 
y(w) — 2\w — 5 crosses the curve y(w) — — ^=t= in two 

points (provided 2A < max{y'(w)} — — 1/\/2), one of the 
points corresponds to the unstable fixed point on the phase 
portraits of Figs.2c-2f, while the other to the stable elliptic 
point at 6 = n in these Figs. As S decreases further, the 
unstable fixed point moves to w = 1. 



II. NONLINEAR TWO-MODE MODELS FOR 
ATOM-MOLECULAR SYSTEMS. 

A. Model equations and its physical origin; 
classical phase portraits 

In BEC-related mean-field models nonlinearity usually 
comes from s-wave interactions (via a scattering length 
entering the nonlinear term of the Gross-Pitaevskii equa- 
tion Q). However, interesting nonlinear models arise in 
atom-molecular systems, where atoms can be converted 
to BEC of molecules. Even neglecting collisions and cor- 
responding s-wave interactions, the nonlinearity comes 
into play from the fact that two atoms are needed to 
form a molecule. 

We consider a Hamiltonian system with the Hamilto- 
nian function 

H = -5(t)w + Xw 2 + (1 - w)y/l + wcos 9, (4) 

where r = ei, e <C 1. Several systems can be described by 
the model dU), in particular coupled atomic and molec- 
ular BEC [2a], and agas of fermionic atoms coupled to 
molecular BEC [H, [23, HH ■ Let us briefly discuss these 
systems. 
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FIG. 4: The Bloch sphere corresponding to ABEC models and the generalized Bloch sphere corresponding to AMBEC models 
(the surfaces u 2 + v 2 = 1 — w 2 on the left and u 2 + v 2 = |(to — l) 2 (f^ + 1) on the right). At large detuning, near w = 1, the 
area within a trajectory on the generalized Bloch sphere is proportional to u 2 + v 2 ~ (1 — w) 2 = T 2 , while on Bloch sphere the 
area is proportional to u 2 + v 2 ~ 2(1 — w) — 2T . Note however that action variable in either case is proportional to 1 — to. 
Action is related to the area on the Hamiltonian phase portraits which is approximately equal to 1 — to for the corresponding 
trajectory, see [27l.[2^|. 



In [351 ] . a system of coupled atomic and molecular 
condensates was considered using generalization of the 
Bloch representation for the two-mode system. Quantum 
Hamiltonian of the system is H = ^-a)a + ^(a^a^b + 
b'aa), where and a are the creation and annihila- 
tion operators of the atomic mode, while W and b are 
the creation and annihilation operators of the molecu- 
lar mode. The two modes are supposed to be coupled 
by means of a near resonant two-photon transition or 
a Feshbach resonance, with a coupling frequency £1 and 
detuning A. Introducing operators L x = \f2 a "J^ " , 

Ly = V2 a ' a ^ aa , L z = 2b V ta , Heisenberg equa- 
tions of motion in the mean-field limit lead to the dy- 
namical system for the rescaled components of general- 
ized Bloch vector: s x — —5s y , s y — — -^(s 2 — l)(3s z + 
1) + Ss x , s z — —\^2s y , where the rescaled detuning is 
S = A/(^/NQ), while s x ,y,z are the expectation values of 
L x ,y,z { s z = 1 corresponds to all-molecule mode; N is the 
number of atoms). Exactly the same dynamical system 
arises in the degenerate model of fermionic atoms cou- 
pled to BEC of diatomic molecules [2(| . Indeed, using a 
similar approach in pol ] an analogous system of equations 
was obtained 

u = S(t)v, 

/2 

v = -5(t)u+^—(w-1)(3w + 1), (5) 
w = V2v, 

where w is the population imbalance between all-atom 



and all-molecule mode, u and v are real and imaginary 
parts of the atom-molecule coherence, 5 is the rescaled 
detuning from the resonance. These equations are equiv- 
alent to the Hamiltonian equations of motion of the 
Hamiltonian system (0| with A = [27|. The vari- 
able 8 canonically conjugated to w is related to the 
old variables as 8 = atan(t>/u). The all-atom mode 
now corresponds to w = 1, while all-molecule mode to 
w = —1. Sweeping through Feshbach resonance from 
fermionic atoms to bose molecules can be described by 
the Hamiltonian (01 with A = and S slowly changing 
from large positive to large negative values. In both sys- 
tems (atom-molecule BEC and degenerate fermionic gas 
coupled to BEC of molecules) mean-field collisional in- 
teractions were neglected so far. The case A ^ in the 
Hamiltonian ^ corresponds to inclusion of the s-wave 
scattering interactions. Recently, in |43j a more general 
quantum Hamiltonian describing the coupling between 
atomic and diatomic-molecular BECs within two-mode 
approximation was considered: 

H = U a Nl + U b N* + U ab N a N b +(i a N a + 

b + tfaa), (6) 

where is the creation operator for an atomic mode 
while b^ creates a molecular mode; parameters Ui 
describe S-wave scattering: atom- atom (U a ), atom- 
molecule (U a b), and molecule-molecule (£4). The pa- 
rameters Hi are external potentials and Q is amplitude 
for the interconvertions of atoms and molecules. N a and 
N b are populations of the atomic and molecular mode, 
correspondingly. In the limit of large N = N a + 2N b , the 
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classical Hamiltonian was obtained: 

H = Xz 2 + 2az + f3 + 2^1^1(1 + z) cos(49/N), (7) 
where 



-(U a /2-U ab /i + U b /8), 



■(U a /2~U b /8 + f, a /2N-fi b /4N), (8) 



is phase difference between the modes, z is difference 
in populations, ft is not difficult to transform the Hamil- 
tonian to the form Q denoting z = — w and intro- 
ducing a new time variable t' = At/N to get rid of the 
4/iV multiplier in the last term of (J7J). The term is not 
important for dynamics. Therefore, the Hamiltonian Q 
describes coupled atomic-molecular BECs (with s-wave 
interactions) in the mean-field limit. Sweeping through 
Feshbach resonance can be modelled now by changing 5 
and keeping A fixed in the Hamiltonian @ . Self-trapping 
phenomenon in the model discussed in [43| allows to pre- 
dict qualitatively new effect, that is non-zero remnant 
fraction in the adiabatic passage through the resonance; 
we do not present detailed quantitative analysis of the 
model in the present paper, but note that it may provide 
an alternative explanation of finite conversion efficiency 
at Feshbach resonance passage within mean-field approx- 
imation. Similar to approach of [26| mentioned above, s- 
wave interactions within molecular BEC can be included 
in the model of fermionic atoms-bose molecules system 
via the same coefficient A ^ in (U). 

Phase portraits with A = (Case I) and different values 
of 8 are given at Fig. [TJ Phase portraits with some con- 
stant A < (Case II) and different values of S are given 
at Figf2] The phase portraits for Case I were analyzed in 
detail in [27|]. The dynamics can also be visualized using 
variables u, v, w of the system ([5]). The latter system pos- 
sesses an integral of motion u 2 +v 2 ~ 1 ) 2 (uj + f ) = 
defining the generalized Bloch sphere (see Fig. 4). The 
important property of the generalized Bloch sphere is 
the singular (conical) point at (0,0,1). As described in 
[221, the points (0, 0, ±1) are represented by the segments 
w = ±1 in the Hamiltonian phase portraits. Neverthe- 
less, it does not mean that all the points of the either 
segment are equivalent. As described in [27], [28j], saddle 
points appear on the segment w — 1 at certain values of 
the parameter 6. This drastically influence dynamics in 
the vicinity of w = 1. Let us briefly recall the description 
of the phase portraits previously given in [27l ] . 

If 5 > v2j there is only one stable elliptic point on 
the phase portrait, at 9 = and w not far from —I [see 
Figure la]. At 6 = v2 a bifurcation takes place, and at 
\[2 > S > the phase portrait looks as shown in Figure 
lc. There are two saddle points at w — 1, cos# = —5/y2 
and a newborn elliptic point at 9 = w. The trajectory 
connecting these two saddles separates rotations and os- 
cillating motions and we call it the separatrix of the 
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FIG. 5: Time evolution of the adiabatic invariant (action) I 
and the improved adiabatic invariant J in the model (|4]) with 
A = 0. 



frozen system (what is most important is that the pe- 
riod of motion along this trajectory is equal to infinity). 
At S — on the phase portrait the segment w = — 1 be- 
longs to the separatrix (Fig. [5]i). At < S < V2 the 
phase portrait looks as shown in Fig. [2^. At 6 = —\[2 
the bifurcation happens, and finally, at large positive val- 
ues of 6, again there is only one elliptic stationary point 
at 9 = 7r, and w close to —1. 

Let us introduce the action variable. Consider a phase 
trajectory on a phase portrait frozen at a certain value 
of 6. If the trajectory is closed, the area S enclosed by it 
is connected with the action I of the system by a simple 
relation S = 2nl . If the trajectory is not closed, we 
define the action as follows. If the area S bounded by 
the trajectory and lines w = 1, = 0, = 2tt is smaller 
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than 2n, we still have S = 2-kI. If S is larger than 27r, we 
put 2itI = 47r — S. Defined in this way, / is a continuous 
function of the coordinates. 

How does the process of Feshbach resonance passage 
happen in terms of the classical portraits of FigQ]7 Sup- 
pose one starts with w(0) — w sa 1, and 6(0) ^> 1 (physi- 
cally, it means that almost all population is in the atomic 
mode, but there is small initial molecular fraction). In 
the phase portrait of the unperturbed system the corre- 
sponding trajectory looks like a straight line (Fig. QJ,). 
The initial action of the system approximately equals to 
1 — wo- For example, assume that the area S* within the 
separatrix loop in Fig. [Tfc (corresponding to 6—6* = 1) 
is equal to S* — 2nI Q — 2n(l — wq). When, as 6 slowly 
decreases, the trajectory on an unperturbed phase por- 
trait corresponding to the exact instantaneous position of 
the phase point {w(t), 0(t)} slowly deforms, but the area 
bounded by it remains approximately constant: action 
is the a ppr oximate adiabatic invariant far from the sep- 
aratrix [l0[. As 6 tends to <5*, the form of the trajectory 
tends to the form of the separatrix loop in Fig. [Tfc. The 
phase point is forced to pass near the saddle point at the 
W = 1 segment many times. Since the area S within the 
separatrix loop slowly grows, approximately at the mo- 
ment t = t* when <5(t*) = 6* separatrix crossing occurs, 
and the phase point changes its regime of motion from ro- 
tational to the oscillatory around the elliptic point inside 
the separatrix loop. Then, it follows this elliptic point 
adiabatically (as no separatrix crossings occur anymore) . 
The elliptic point reaches w = — 1 at large positive 6. 
The value of the population imbalance tends to some 
final value w = wj. The action variable at large 6 is ap- 
proximately equal to 1 + w. We see that in the adiabatic 
limit the sign of the population imbalance is reversed, 
wo = —Wf. Nonadiabatic correction to this result arise 
due to the separatrix crossing and is discussed in detail 
in the next paragraph. 

In the Case II the phase portraits have richer structure 
(Fig. [5]). With A ^ 0, another saddle point can appear 
at 9 = 7r provided A < A c = —5572-. The appearance of 
this saddle point can be understood from the graphical 
solution of the equation (see also 
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As 6 is decreased, the line y(w) = 2Xw — 6 goes up and 
crosses the curve determined by the r.h.s of ©. Two 
points of intersection represent the saddle point (which 
moves to w — 1 as <5 is decreased further) and the elliptic 
fixed point which moves to w = — 1. As the saddle point 
reaches the w = 1 segment, another bifurcation occurs 
and the saddle point "splits" into the two saddle points 
(similar to those in Fig[T]), that move apart from 9 = n 
along the segment w = 1 and disappear at 9 = 0. 

In Section lib the dynamical change in the action in 
the case A = is considered in detail, while Section lie 
briefly discusses the case A 7^ (geometric jump). 



B. Case I: negligible mean-field interactions, A = 0. 
Dynamical change in the action at the separatrix 
crossing. 

Consider in a greater detail the passage through the 
separatrix in Fig. [T]described in the previous subsection. 
At large positive 6, 1 — w is proportional to classical ac- 
tion, while at large negative 6 action is proportional to 
1 + w (see also Fig. [4j. In the adiabatic limit, w reverses 
its sign due to passage through the resonance: the final 
and initial values of w are related as uif = —uii n . Cal- 
culating c hang e in the action due to separatrix crossing 
(Refs. [27, 28]), one obtains the nonadiabatic correction 
to this adiabatic result. It scales linearly with e if initial 
population imbalance slightly deviates from 1 (i.e., initial 
molecular fraction is not very small). 

As the trajectory nears the separatrix due to slow 
change (of order e)in the parameter, the action undergoes 
oscillations of order of e. Each oscillation corresponds to 
one period of motion of the corresponding trajectory in 
the unperturbed system. In the vicinity of separatrix, 
the period of motion grows logarithmically with energy 
difference h between energy level of the unperturbed tra- 
jectory and the energy on the separatrix (so as h tends to 
0, the period of motion tends to infinity). As a result, the 
"slow" change of the parameter becomes "fast" as com- 
pared to the period of motion: breakdown of adiabaticity 
happens; oscillations of the adiabatic invariant grow and 
at the crossing its value undergoes a quasi-random jump 

(Fig. 0. 

According to the general theory, it is not enough to 
consider dynamics of the action variable. One introduces 
the improved adiabatic invariant J = I + ef(w, 9, t) (see 
the Appendix for brief description of adiabatic and im- 
proved adiabatic approximations and the general formula 
for J ). The improved adiabatic invariant is conserved 
with better accuracy: far from the separatrix, it under- 
goes very small oscillations of order e 2 . At the separatrix 
crossing, it undergoes jump of order e. 
We illustrate this behavior in Fig. [5) Figs. [5^,,b give dy- 
namics of the action (adiabatic invariant) /. It is clearly 
seen that before and after separatrix crossing it oscillates 
around different mean values, but the jump in action is 
of the same order as its oscillations close to the separa- 
trix. Fig. presents time evolution of the improved 
adiabatic invariant. The jump in J is much more pro- 
nounced (although it is possible to express the improved 
adiabatic invariant in the elliptic functions, we choose to 
calculate it numerically according to the definition given 
in the Appendix A). 

Now, at large \6\ not only the action I coincides with 
value of 1 — \w\, but also the improved adiabatic invariant 
J coincides with /. Therefore, calculating change in the 
improved adiabatic invariant J, we obtain change in the 
action and change in the value of 1 — \w\ due to the 
resonance passage. For the case of small initial action J, 
the change in action was calculated in [27j according to 
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FIG. 6: Scattering at the separatrix crossing, a) Bunch of trajectories with various (but close) initial conditions undergoing 
jump of the improved adiabatic invariant at separatrix crossing. Trajectories are mixed due to the jumps, b) e— dependence 
of magnitude of jump of the improved adiabatic invariant. For every value of e, we calculated a bunch of 80 trajectories from 
8 = 10 to 8 = 0. Initial values of w were chosen to be equidistantly distributed in the interval [0.96, 0.96 + 1.5e]. The theory 
predicts quasi-random jump of the improved adiabatic invariant, which magnitude scales linearly with e. We calculate mean 
value a of squared change in the improved adiabatic invariant, which turns out to scale perfectly linearly with e (accordingly, 
dispersion a 2 scales linearly with e 2 ) c) The same as in Fig.b, but with smaller values of e and initial values of action. We 
calculated slope of the line cr(e) taking into account the four points with the smallest values of e, and get the value k w 1.1614, 
which is in good correspondence with theoretical prediction of -y/4/3 fa 1.15; for larger values of e, correspondence worsens: 
k ~ 1.015 when taking into account all points, d) High sensitivity of the jump of the adiabatic invariant on initial conditions. 
Calculations for e = 0.0004 are presented. Initial values of w for 100 trajectories were uniformly distributed in the tiny interval 
(wo,wo + 1.5e). Change in the improved adiabatic invariant was calculated (AJ = J(S = 0) — J (8 = 10)). It is seen that tiny 
change in the initial conditions results in large variance of the jump of the action. Trajectories arrive at the separatrix with 
different values of the pseudo-phase £. Maxima in the Figure correspond to £ = and £ = 1. The formula for the jump of the 
adiabatic invariant predicts high increase in the value of the jump when ~ (tt£) nears 0. In the very vicinity of £ = 0, 1 the 
formula is not working (the predicted jump diverges while the calculated jump is finite), however measure of the exceptional 
initial conditions leading to £ = 0, 1 is very small [10| . 



the general method of The formula is 



2ttAJ 



•> e€> * : hi(2 sin 7Tg), 



(10) 



where is rate of change of the area within the sepa- 
ratrix loop: = (note that the rate do not depend 
on e); £ is the pseudo-phase: £ = \ho/eQ\, where ho 



is the value of the energy at the last crossing the ver- 
tex bisecting the angle between incoming and outgoing 
separatrices of the saddle point C outside the separatrix 
loop (see FigOJ). Similar calculations were done in [2fj| . 
The formula can be further simplified by expressing 
via || = 5' '. Indeed, the area within the separatrix loop 

- arccos ( ) . We are inter- 

^ Vl+w ) 



is S(5) = 2 J 5 1 2 _ 1 dw 
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ested in derivative of S(8) over S. Differentiating the 
above integral over parameter <5, one obtains: 

at 

S\6) = 4^2 -<5 2 , and 6 = S'(5)— = S'{6)6'. 

OT 

(11) 

Therefore, the formula (fTO]) is simplified to 

AJ = -— In(2sin7r0. (12) 

We carefully checked numerically behaviour of jumps 
in action predicted by formula (fT2"j) . see Fig.©. Fig[6ji 
demonstrates scattering at the separatrix crossing: 
bunch of trajectories with various (but close) initial con- 
ditions undergoing jumps of the improved adiabatic in- 
variant at separatrix crossing. Figures [6Jd, c demonstrate 
e— dependence of jumps of the improved adiabatic invari- 
ant. For several values of e, a bunch of 80 trajectories 
with close initial conditions was calculated, and disper- 
sion of actions due to separatrix crossing was calculated, 
which scales linearly with e 2 (i.e., a 2 — Ke 2 ). Note that 
from the formula (fT2")) it is possible to determine not only 
the linear power-law, but also the corresponding coeffi- 
cient of proportionality K . The theory predicts uniform 
distribution of £, therefore dispersion of jump in the ac- 
tion can be calculated as 

a 2 = 16e 2 (5') 2 ir- 2 [ In 2 2 sin tt^ = (13) 
Jo 3 

For numerical calculations, we used linear sweeping of 
6, therefore the predicted dispersion is a 2 = 4e 2 /3. Pre- 
dicted coefficient 4/3 can be compared with the slope in 
Figs[6j3,c . For relatively large e (Fig. [6Jd), correspon- 
dence is not very good, but when we decreased the value 
of e, we obtained K 1.348 which is in good correspon- 
dence with theoretically predicted K = 4/3 w 1.333. 
We reveal also high sensitivity of the jump of the adi- 
abatic invariant on initial conditions (Fig[SJi), which is 
the cause of uniform distribution of £ [H| . We therefore 
checked almost all qualitative and quantitative aspects 
of destruction of adiabatic invariance at separatrix cross- 
ings in that model. Let us finally mention the main steps 
in obtaining the formula : 

1. Linearization around the saddle point in the frozen 
system and derivation of a approximate formula for 
the period of motion T along the trajectory with 
energy h. The period depends logarithmically on 
h and is inversely proportional to the square root 
from the Hessian of the Hamiltonian in the saddle 
point (determinant of the matrix of second deriva- 
tives) . 

2. Obtaining the action variable / from the period T 
using the formula T = 2%dl/dh. 

3. Calculation the function / at a point of the vertex 
bisecting the angle between incoming and outgoing 



separatrices of the saddle point (Fig. lc). It is 
proportional to O (for details, see [lOj). 

4. "Slicing" the exact trajectory on parts (correspond- 
ing to "turns" in the unperturbed system) by the 
bisecting vertex and constructing a map r„, J n — > 
T n +i i Jn+x using the analysis described above (tq 
is the moment of last crossing of the vertex before 
the separatrix crossing, t_i is a previous moment 
of crossing the vertex, etc. J n is value of the im- 
proved adiabatic invariant at r„). Summation of 
changes of adiabatic invariant at each turn leads to 
the formula (fT2"|) . 

See Refs. [H [H| for further details. 



C. Case II: Condensates with interactions, A / 0. 
Analog of nonzero adiabatic tunnelling. 

Let us briefly consider the model with A < A c = — 7^72 • 
Separatrix crossing happens via another scenario here 
(according to the motion of the fixed points described 
in Section Ha) . We give only qualitative discussion of a 
possible new phenomenon. We plot the phase portraits 
at different 5 and fixed A in Fig. @. Now, as S is de- 
creased, three domains Gi 2,3 appear in the phase por- 
trait at certain r = as a result of the first bifurcation. 
Shortly after the bifurcation (see Fig. [2):) the separatrix 
consists of the two "loops": the upper G2 (adjacent to 
w = 1 line), which area 52 (r) decreases from 52 (r*) to 
zero as the unstable fixed point goes towards w = 1, and 
the bottom G3, whose area 53 (r) increases from zero; G\ 
is the "outer" domain adjacent to w = — 1. 

In case initial action Iq of a phase point is sufficiently 
small (2ttIq < 52(r*)), the phase point resides in the 
G2 domain when the separatrix emerges (without any 
separatrix crossing, see Fig[2fc). In case 2ttIo is larger 
then the area 52 of the domain G2 at the moment of 
separatrix creation, the phase point occupy G\ at this 
moment. Consider the former case, i.e. small initial ac- 
tion. As 5 evolves, 52 decreases, while 53 grows. When 
52 (t) becomes equal 2ttI , separatrix crossing occurs and 
the phase point is expelled to G\ domain and then al- 
most immediately to G3 domain (say, in the Fig. [5f). It 
is easy to see that the phase point acquires large action 
due to geometric jump in the action when entering G3, 
so in the end w will deviate from the all-molecule mode 
w = — 1 considerably (the geometric jump is equal to 
(53(1-**) — 52(t**))/27t, where r*» is the moment of the 
separatrix crossing). This is in some sense analogous to 
the nonzero A T discussed in [24l.[25j and in Section III of 
the present paper. One might try to explain the sizable 
remnant fraction after the adiabatic Feshbach resonance 
passage as the geometric jump in the action due to the 
self-trapping effect of s-wave interactions. This, however, 
requires further investigation: while calculation of the ge- 
ometric jump in action is trivial, dynamical jump is not 
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so easy to calculate in this geometry. So far, we just sug- 
gest a possible new phenomenon in the model, detailed 
discussion will be given elsewhere. 



III. NONLINEAR TWO-MODE MODEL FOR 
TWO COUPLED BEC. 

A. Model equations and its physical origin; phase 
portraits 

We consider the Hamiltonian ("nonlinear 2-mode 
ABEC model") 

Xw 2 



H = -6w 



(14) 



There are many systems in BEC physics that are de- 
scribed in the mean-field limit by the Hamiltonian (|14p . 
It has been used to model two coupled BECs (BEC in 
a symmetric double well in case 6 = 0) ([23|). The 
model with 8 ^ is equivalent to nonlinear Landau- 
Zener model, which appears, in particular, in studying 
BEC acceleration in optical lattices [13, HE] 

Theory of nonlinear Landau-Zener tunnelling was sug- 
gested in [2^. [2511 . However, only the case of zero initial 
action was considered. In particular, it is said in [25j 
that adiabaticity is broken when "fixed points collide". 
In case the initial action is not zero, adiabaticity is broken 
before that: it is broken when separatrix crossing occurs. 
Therefore, it is necessary to involve theory of separatrix 
crossings in consideration of these models. 

It is worth to mention that for BEC in a symmetric 
double- well, there exist also improved 2-mode model [38[ , 
where the term cos 20 is added: 

z 2 1 

H = A— - B a/1 - z 2 cos 4> + -C(l - z 2 ) cos 20, (15) 

where parameters A, B,C are determined by overlap in- 
tegrals of the mode functions. Usually, the cos 20 term 
is small and can be omitted. Then, the improved model 
Hamiltonian can be reduced to (Q3J) with 8 = (still, 
coefficients are determined more accurately in the im- 
proved model). The original model is derived for the case 
of constant parameters. One may wonder if it is work- 
ing in a time-dependent situation. It is not difficult to 
demonstrate that for slowly changing parameters one can 
use the same model, with parameters of the Hamiltonian 
slowly changing in accordance with the "instantaneous" 
model. For simplicity, let us demonstrate this using the 
improved 2-mode model [38| as an example. The order 
parameter in a two-mode approximation is 



ip(x,t) 



$+(x) ± $(X) 



(16) 



$i, 2 (x) 



V2 



where $± satisfy the stationary CP equation 
ld 2 $± 



The variables of the classical Hamiltonian are defined 



z(t) = \Mt)\ 2 - \Mt)\ 2 , <P(t) = arg^W - arg^i(t) 

(18) 

Substituting (TiUl) , Q17p into the time-dependent GP equa- 
tion, one gets [38j 



z±(M*)±Mt))\P± - gN\$ ± \ 2 }<S>± + 



(19) 



$ 2 ± $ T Q ± ] 



where P±,Q± are functions of ipi,ip2 (see [38|]). From 
these equations, one get the equations of motion for 
ipi,tp2 (Eqs. 13 from [381]): 



= (F + A\^ 



A7 



A/3 £7 2 

( — 2" + t'^' 



(20) 



#fc*± 



2 dx 2 



U ext $ ±+5 |$ ± | 2 $ ± (17) 



which can be rewritten as Hamiltonian equation 
of motion of the corresponding classical pendulum 
(F, A, C, A7, A/3 are functions of mode overlap integrals 
and energies (3±). Considering time- varying parameters, 
we introduce instantaneous mode functions <f>±(x,t). If 
we keep two-mode expansion of the order parameter, 
when it is not difficult to show that additional terms 
coming from time-dependence of the mode functions 
(f $ + ^-dr, J <& + ^-dr, etc ) are strictly zero due to 
symmetry and normalization conditions. Complications 
can arise only from excitation of other modes (if we would 
allow, say, four-mode expansion). However, we do not 
consider this question here. Even in the two-mode ap- 
proximation nonadiabatic dynamics is nontrivial, and it 
comes purely from nonadiabatic behaviour of classical ac- 
tion. Phase portraits of the model (fT4"|) with 6 = are 
given in Fig. [7] We are interested only in the supercrit- 
ical case here. Separatrix crossings and corresponding 
changes in the action are discussed in Section lib. The 
case 8 (NLZ model) is discussed in Section lie, where 
we present a new phenomenon: separated AT. 



B. Case I: symmetric double-well, 5 = 0. 

We suppose initially the system is in the oscillating 
regime of complete tunnelling osciallations (domain G3), 
and then due to slow change of parameters is switched 
into self-trapped regime. Two different probabilistic phe- 
nomena take place at the crossing: quasi-random jump 
in the action and the probabilistic capture. 

Consider the probabilistic capture: there are two do- 
mains G12 for the self-trapped regime in the phase por- 
traits: in the first (upper) w > 0, in the second (bottom) 
w < 0. In which of these two domains the phase point 
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FIG. 7: Phase portraits of the 2-mode ABEC Hamiltonian with 5 = 0. From left to right (a,b,c,d): A = 20,2.4, 1.2,0.8. As A 
decreases, separatrix loop grows until A = 2 where it changes its configuration, and at A = 1 it disappears. On the other hand, 
by increasing A it is possible to switch from regime of complete oscillations (domain 3) to the self-trapped regime (domains 1 
or 2). The unstable fixed point do not move: it is either at (0,0) or absent. 




-1 0.5 0.5 1 

w 



FIG. 8: Graphical solution of the equation —S + Ait; = 
w/y/1 — w 2 which gives fixed points at 8 — ir. As S decreases, 
the line goes up, and three fixed points can appear from a sin- 
gle one at certain window of value of 5 provided A > 1. The 
star denotes the unstable fixed point which after the birth 
goes down and collides with the stable fixed point. See corre- 
sponding phase portraits in the next Figure. 



[10j. As a result, the formula for the jump of the action 
is simplier than expressions for the action itself. Suppose 
A > 2 so the phase portrait looks like in Fig. [7Jd and 
we start from the regime of complete oscillations. Slowly 
changing A, we can switch to the self-trapped regime. 
The expression for the area of the separatrix loop is easy 
to calculate: 



S(t)/4 = b + arcsin 6, b = (21) 

A 

The Hessian of the Hamiltonian in the unstable fixed 
point can be calculated as D(t) — — (A — 1). 
Let us define 

d(r) = (22) 

We calculated jump of the action according to the gen- 
eral method as 



will be trapped (in other words, in the left or the right 
well)? The trapping in either of the domains is also very 
sensitive to initial conditions; in the limit of small e the 
trapping is a probabilistic event. For the symmetric case, 
the probability to be trapped in ether well is exactly 1/2. 
However, for the asymmetric well the answer is not so 
straightforward. It is determined by some integrals over 
separatrix at the moment of switching (general theory 
exists, see (Toj). 

As for the first phenomenon (jump in the action), at 
the moment of switching, destruction of adiabaticity hap- 
pens in the sense that the adiabatic invariant undergoes 
a relatively large jump of order of i/e (very similar to 
that discussed in the Section II). If we then slowly bring 
the parameters back to the initial values, the adiabatic 
invariant will be different. 

The formulas for the action-angle variables are cum- 
bersome. In fact, to calculate change in the action, it is 
not necessary to have formulas for the action-angle vari- 
ables. The jump is determined by local properties of the 
Hamiltonian near the separatrix: the area of the sepa- 
ratrix loop and the Hessian of the unstable fixed point 



1 4A' 

A J = ed,6* ln(2 sin«)) = e— ■ ln(2 sin«)), 

27T n\ z 

(23) 

where £ is the pseudophase corresponding to the first 
crossing of line 9 = n in the G\ y 2 domains, d» is value of 
d at the moment of crossing this line, values of A and A' 
are also taken at this moment. 

We checked this formula numerically. A set of 100 
phase points with initial conditions being distributed in 
a small (of order e) interval far from the separatrix were 
chosen. Then, the bunch of trajectories in the system 
with slowly changing parameter was calculated. For each 
trajectory, values of £ and A J (change in the improved 
adiabatic invariant) were determined. From numerically 
determined £, theoretical prediction for change in the ac- 
tion AJ was calculated and compared with numerically 
determined AJ. Results are in the Fig. [101 correspon- 
dence between numerical results and analytical predic- 
tion is perfect. In the same calculations, mechanism of 
quasi-random division of phase flow was verified: exactly 
one half of the phase points from the considered set were 
captured in the upper domain G\, and the other half 
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FIG. 9: Nonlinear Landau-Zener tunneling: phase portraits of the 2-mode ABEC Hamiltonian at different values of 5. From 
top left to bottom right: 5 = 20, 3, 1.8, 1.2, 0, -1.2, -1.8, -3, -20; A=const=4. 



were trapped in the lower domain Gi. This is a purely 
classical phenomenon, the sound example of probabilistic 
phenomena in dynamical systems ( [T3, !H| ) . 



C. Case II: asymmetric double- well and nonlinear 
Landau-Zener model, 8 7^ 0. Separated adiabatic 
tunnelling. 

Consider sweeping value of 8 from larg e positive to 
large negative values in Fig{9] (see also [2J, [25|). In case 
A < 1, only two fixed points exist at 9 — 0, ir (P2, Pi cor- 
respondingly). As 8 changes from 8 = — 00 to 8 = +00, 
Pi (corresponding to the lower "eigenstate") moves along 
the line 9 = ir from the bottom (w = —1) to the top 
(w = 1), the other point Pi (corresponding to the up- 
per "eigenstate") moves from the top to the bottom. In 
case A > 1, two more fixed points appear in the window 



-8 C < 8 < 8 C , 8 C = (A 2 / 3 - I) 3 / 2 . We concentrate on 
this, "above-critical" case. The new points lie on the line 
9 = 7r, one being elliptic (P3) and the other hyperbolic 
(P4). Again, it is convenient to use graphical solution 
(Fig. |HJ) to visualize appearance and disappearance of 
the fixed points. It is stated in the (25[, that collision be- 
tween Pi and P3 leads to nonzero AT from the lower level 
to the upper level, and tunnelling probability in the adi- 
abatic limit is obtained by calculating phase space area 
below the "homoclinic trajectory" ( which is the limiting 
case of the separatrix with S3=0), i.e. as geometric jump 
in the action. In the zeroth order approximation, this 
approach is correct (if the initial action is zero or very 
small). It is very important that we can adopt general 
theory of separatrix crossings to the case of this model 
with nonzero initial action (corresponding to initially ex- 
cited system). In this case destruction of adiabaticity 
happens before collision of the phase points. Initial tra- 
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FIG. 10: Jump in the improved adiabatic invariant in dependence of the pseudo-phase £. Filled squares: numerical results; 
dashed line: analytical predictions according to Eq. 1231 We slowly changed A according to the law X — \ a — \b cos et, with 
e = 0.001, X a = 15, Xb — 10. We took a set of 100 phase points with different but very close initial conditions: Wi = 0, 9i are 
distributed along an interval of order of e at the time r = et=0. We propagate the bunch of trajectories until the time r = n 
(so all the points changed its regime of motion from complete oscillations to the self-trapped mode). For each point, value of £ 
and change in the improved adiabatic invariant AJ was determined numerically, then the analytical prediction for the change 
in the improved adiabatic invariant AJ(£) was calculated according to Eq |23l The numerical and analytical results shown in 
the Fig. (a) are almost indiscernible, in (b) enlarged part of the same plot is presented, where small deviations are seen. It is 
also important to mention, that from 100 phase points exactly 50 were trapped in the upper domain Gi, and 50 in the lower 
G 2 . 



jectory is almost a straight line, so the initial action is 
equal to w + 1 in case we start close to w — — 1, or 1 — w 
in case we start close to w = 1. Consider the former 
case. Let the initial action I (i.e., value of w + 1 in 
Fig. [5^) be equal to area of the separatrix loop in Fig. 
Eg. The phase point is oscillating around slowly moving 
P\ point until the area of the separatrix loop S±(t) be- 
comes equal to 2ttIq at some moment t = r*. Where, 



separatrix crossing occurs. Action undergoes geometric 
jump (which is simply [S^t*) — Si(t*)]/2tt. This geo- 
metric jump is the analog of AT probability discussed in 
[24l [25j for the case of zero initial action. The geometric 
jump is accompanied by the dynamical jump similar to 
that discussed in Section II and Section Illb. The dy- 
namical jump is small (of order of e) as compared to the 
geometric jump. But conceptually it is very important: 
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only dynamical jump leads to destruction of adiabatic in- 
variance in the model. Indeed, if we reverse change in 6, 
the phase point will return to its initial domain and the 
geometric jump will be completely cancelled. However, 
dynamical jumps will not be cancelled, and at multiple 
separatrix crossings they lead to slow chaotization (see, 
for example, [l||). Formulas for the dynamical jumps in 
the asymmetric case are more complicated, as there are 
terms of order e and e lne. However, the probabilistic 
capture in this case is very much different. Consider the 
phase portraits in Figs. [5f,g. Suppose that not only 5, 
but also A is changing. At the moment of crossing, the 
area S3 is diminishing, while the areas can behave 
differently depending on evolution of parameters. Sup- 
pose both 5*1,2 are increasing: 61,2 > 0, 83 < 0. De- 
note as Z12 the parts of the separatrix below and above 
the saddle point, correspondingly. There is phase flow 
across I2 from the domain G2 to G\, and across l\ from 
G3 to G2. The latter flow is divided quasi-randomly be- 
tween G2 and G\ : the phase point leaving G3 can remain 
in G2 or be expelled to G2. This is "determined" dur- 
ing the first turn around the separatrix. After that, the 
particles are trapped either in G\ or G2. Probability for 
either event can be calculated as integrals over the sepa- 
ratrix parts hfl ( (Iol|): 

Pi = V 2 = (24) 

r/c ^ / / , fdH dH s \ 

Here integrals are taken along the unperturbed trajecto- 
ries at the moment of separatrix crossing (or last cross- 
ing the line 9 = ir before the separatrix crossing), H s 
is the (time-dependent) value of the Hamiltonian H in 
the unstable fixed point, H denote the Hamiltonian H 
normalized in such a way as to make value of the new 
Hamiltonian in the unstable fixed point to be zero. It is 
possible to calculate all the integrals analytically, see the 
Appendix B. We present numerical example in Fig. [11] 
A set of N = 100 trajectories was considered with initial 
conditions distributed in a tiny interval of w, and with 
9(0) = (so initial actions were distributed in a tiny in- 
terval of order e: Ik = h + kSI, NSI ~ e, k = 1, .., N; 
alternatively, one can consider a set of phase point with 
equal initial actions, but with distribution of phase along 
2ir interval). Both 5 and A were changed; so after the 
separatrix crossing a phase point can be trapped either 
in Gi or G2. From the set of 100 points, 87 were trapped 
in Gi, while 13 were trapped in G2. The difference be- 
tween the final actions of these two subsets is approxi- 
mately Iq, the initial action of points in the bunch. The 
probability of 87% is in good correspondence with the 
theoretical prediction, which gives P 2 = 86, 998 for the 
probability of capture into the domain G2. Possible ex- 
perimental realization of this new phenomenon is again 
BEC acceleration in optical lattices, but with simultane- 
ous modulation of the lattice potential depth. 



IV. CONCLUSION 

We discussed destruction of adiabatic invariance in sev- 
eral nonlinear models related to BEC physics. We espe- 
cially concentrated on the cases that have not been con- 
sidered in the corresponding papers on BEC dynamics 
yet: that is, when the initial action is not zero. 

We found that the general theory of adiabatic separa- 
trix crossings works very well in the considered models. 
Two aspects of destruction of adiabatic invariance were 
considered: quasi-random jumps in the approximate adi- 
abatic invariants and quasi-random captures in different 
domains of motion at separatrix crossings. 

We discussed quasi-random jumps in the approxi- 
mate adiabatic invariants in the models describing Fesh- 
bach resonance passage in coupled atom-molecule BECs, 
BEC tunnelling oscillations in a double well, and nonlin- 
ear Landau-Zener tunnelling. Comparingwith previous 
analysis of the abovementioned models [25|, [26| , the key 
feature of our approach should be emphasized: the sys- 
tem is linearized near the hyperbolic fixed point, not near 
elliptic fixed points of the unperturbed system. 

Another important class of phenomena considered here 
is probabilistic captures into different domains of motion. 
They were discussed for the case of BEC tunnelling oscil- 
lations in a (symmetric or asymmetric) double-well and 
the NLZ model with time-dependence of the nonlinearity 
A. Separated AT was discovered in the latter case. We 
suppose it can have experimental applications in BEC 
manipulations with optical lattices. The conceptual phe- 
nomenon of probabilistic capture was firstly discovered 
in celestial mechanics (while studying resonance phenom- 
ena in Solar system) . It is interesting to draw an analogy 
between intricate phenomena of celestial dynamics and 
phenomena happening in many-body quantum systems. 
Conceptual phenomena related to the classical adiabatic 
theory (which includes both adiabatic invariants and the 
adiabatic (geometric) phases) have recently become im- 
portant trend of research in the highly interdisciplinary 
BEC physics field (see (HHlMllof). We hope the com- 
prehensive analysis presented in this paper adds impor- 
tant contribution to understanding nonlinear dynamics 
of Bose-Einstein condensates. 
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FIG. 11: (Color online) Separated AT. We took a bunch of 100 trajectories with initial actions distributed in a tiny interval 
of order e (in the Figure: the bottom curve consists of 100 initial trajectories; the two upper curves consist of 87 and 13 
trajectories, correspondingly, and represent "snapshot" of the trajectories after the end of sweeping parameters). As parameters 
are changing, the separatrix (see Figs. 9d-f) moves down, and crosses the bunch. Due to quasi-random division of phase flow 
described in the text, some of the points were captured to the Gi domain, while the other to the G2 domain. As a result, phase 
points undergo different geometric change in the action. After the capture, actions are conserved. Therefore, a phase point 
can acquire two different values of the adiabatic invariant. The difference between the values corresponds to the area of the 
domain G3 at the moment of separatrix crossing, i.e. it is approximately equal to the initial action. The question is how the 
initial bunch is divided, what is the probability for a phase point to come into either of the two upper bunches. From the set 
of 100 points, 87 were trapped in the upper bunch, while 13 in the bottom. This numerical result is in very good accordance 
with the theoretical prediction for the probabilities (|24I37[) . which gives Vi = 86.998 (see the Appendix B). In case of nonzero 
AT considered in [2J, [2^|, initial bunch would lie near w — — 1 line, and there would be only one "final" bunch. Here, there 
are separated bunches, which suggested us to introduce the terminology "separated AT". 
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VI. APPENDIX 

A. Adiabatic and improved adiabatic 
appr oximat ions 

To consider change in the action during a separatrix 
crossing, it is necessary to introduce improved adiabatic 
invariant J in addition to the ordinary action variable /. 
Improved adiabatic approximation is discussed in [l(|- 

Let I = I(w,8,t), <j> — cf)(w,6,T) mod 2ir be 
the action-angle variables of the unperturbed (r=const) 
problem. The "action" I(w,9,t) multiplied by 27T is the 
area inside the unperturbed trajectory, passing through 
the point (w, 9) (provided the trajectory is closed; other- 
wise the area of a domain bounded by the trajectory and 
lines 8 — 0, 2ir is calculated). The "angle" is a coordi- 
nate on the same unperturbed trajectory. It is measured 
from some curve transversal to the unperturbed trajecto- 
ries. The change (w, 9) — > (/, <f) is canonical (and can be 
done using a generating function which depends on t). In 
the exact system (with f = e 7^ 0) the variables / and <f> 
are canonically conjugated variables of the Hamiltonian 



where Hq(I,t) is the initial Hamiltonian E(w,9,t) ex- 
pressed in new variables, while the perturbation Hi 
comes from the time derivative of the generating func- 
tion. In case the angle (f> is measured from some straight 
line <f> = const, one has the formula [Io| 

1 [* (BE /dE\\ ,, 8H 

where the brackets < .. > denote averaging over the "an- 
gle" 0. _ 

Consider a phase point of the exact system with the 
initial conditions I = Iq, <j) = <^o- The adiabatic approxi- 
mation is obtained by omitting the last term in (|25|) and 
gives 

1 r* 

I = I Q , 0=0 O + - / UJ (I,T)dT (27) 

e Jo 

. Improved adiabatic approximation is introduced in the 
following way. One makes another canonical change of 
variables (I,4>) — > (J, i/j). The change is 0(e)— close to 
the identity and in the new variables the Hamiltonian 
has the form 

H = H (J,t) +eH 1 (J,r) +e 2 H 2 (J,Lu,T,e), (28) 

H, =< « >= -i f (l - f) f * (29) 
ujo Jo V 2 27r / dr 



H = H (I,t)+cH 1 (I,(P,t), 



(25) 



The improved action variable can be defined as 
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J = J{w,9,t) +I + eu, 

If 



u = u(w,6,r) 



dE 



dt, 



(30) 
(31) 



where the integral is taken along the unperturbed trajec- 
tory passing the point (w, 9), T = ^ is the period of the 
trajectory, and the time t is measured starting from the 
point (w, 9). Determined in this way, < u >= 0. The im- 
proved adiabatic approximation is obtained by omitting 
the last term in ([29]) and gives 



J= J , V = "00 + - / (w (J,r) +GJ 1 (J,T))dr, 



dFi 
~dJ' 



(32) 



B. Probabilities of captures during separated AT 

We change both S and A linearly in time: 5 ~ Sq — et, 
A = A — net, k = 1.5; A = 25, Sq = 8. We consider 
a bunch of N = 100 trajectories with initial conditions 
Wk = wq + O.Q2ek, 9k — (wq — —0.8) which imply 
distribution of initial actions in a tiny interval of order 
e. Alternatively, one can consider initial conditions with 
the same initial action, but with distribution along the 
angle variable <f>. In any case, from N trajectories, ap- 
proximately Vi N will be captured in domain G2, and 
V\N in domain G\. As a result, after sweeping value of 
5 to —00, one obtains two bunches of trajectories each 
closely distributed along two different values of action. 
This is a new phenomenon in the context of nonlinear 
Landau-Zener tunnelling. 

At the moment of separatrix crossing, phase portrait 
looks like shown in Fig. [9j. Phase flow from the domain 
G3 is divided between G\ and It is possible to cal- 
culate analytically the probabilities of captures in either 
domain. The separatrix crosses the line 9 = at points 
w = i"a,6, Wa < wi, and the line 9 = it at w — w s (the 
unstable fixed point). These three magnitudes (w ai b,s) 
are the roots of the equation 



(w) 2 = l-w 2 -(h s + 5*w - yw 2 ) 2 = 0, (33) 

where h s is the energy on the separatrix at the moment 
of crossing, and 5*, A* are values of the parameters at 



this moment (w = w a is the doubly degenerate root). In 
other words, 



A2 

w = ±\j — -±(w - w a ){w - w b ){w - w s ) 2 (34) 



Probabilities of capture in either domain are given by 



V2 = -f, V x 



1 P f)M \ f 

w — w« X' f w " , w 2 — w 2 



5' / dw 



+ 



dw- 



2 Jw a ,„ "' 



where lower limits of integration for 1\ , J 2 are w a and 
Wb correspondingly. For value of w one uses the Eq. [34] 
which makes the integrands in Eqs. [35] simple, and one 
gets 



-r 



I 1 = 



-r 



-2w s + w a + w b 



w b - w a 



vr/2, 



y/-{w s ~ W a ){w s - Wb) + (w s + (w a + Wb)/2) l(, 
—2w s + W a + W b 



W b - W a 



tt/2, 



(36) 



- y/-{w s - w a )(w s - w b ) + (w s + (w a + W b )/2)L 



Therefore, 



h = -S'{-a - tt/2) + \[-Q s + W s (-a - tt/2)] 
h 



a 



~8'{a - tt/2) + f [Q s + W s (a- tt/2)} 
-2w s +w a + w b 



W b - W a 



(37) 



y/-{w s - W a )(w s - Wb), W s — w s + (w a + w b )/2 



In the numerical example presented in Fig. 1111 5' = 
— 1, A' = — k = —1.5; at the separatrix crossing A* = 
8.3863369, 5* = -3.0757753, h s = 0.3553544. It gives 
w a w -0.9239628, w b » 0.30155167, w s « -0.4223149. 
The formula (|37|) gives V2 ~ 86.998, which perfectly cor- 
responds to the numerical result (87%). 
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